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Vesicles and many biological membranes are made of two monolayers of lipid molecules and form 
closed lipid bilayers. The dynamical behaviour of vesicles is very complex and a variety of forms 
and shapes appear. Lipid bilayers can be considered as a surface fluid and hence the governing 
equations for the evolution include the surface (Navier-)Stokes equations, which in particular take 
the membrane viscosity into account. The evolution is driven by forces stemming from the cur¬ 
vature elasticity of the membrane. In addition, the surface fluid equations are coupled to bulk 
(Navier-)Stokes equations. 

We introduce a parametric finite element method to solve this complex free boundary problem, and 
present the first three dimensional numerical computations based on the full (Navier-) Stokes system 
for several different scenarios. For example, the effects of the membrane viscosity, spontaneous 
curvature and area difference elasticity (ADE) are studied. In particular, it turns out, that even in 
the case of no viscosity contrast between the bulk fluids, the tank treading to tumbling transition 
can be obtained by increasing the membrane viscosity. Besides the classical tank treading and 
tumbling motions, another mode (called the transition mode in this paper, but originally called 
the vacillating-breathing mode and subsequently also called trembling, transition and swinging 
mode) separating these classical modes appears and is studied by us numerically. We also study 
how features of equilibrium shapes in the ADE and spontaneous curvature models, like budding 
behaviour or starfish forms, behave in a shear flow. 

PACS numbers: 87.16.dm, 87.16.ad, 87.16.dj 


I. INTRODUCTION 

Lipid membranes consist of a bilayer of molecules, 
which have a hydrophilic head and two hydrophobic 
chains. These bilayers typically spontaneously form 
closed bag-like structures, which are called vesicles. It is 
observed that vesicles can attain a huge variety of shapes 
and some of them are similar to the biconcave shape of 
red blood cells. Since membranes play a fundamental role 
in many living systems, the study of vesicles is a very ac¬ 
tive research field in different scientific disciplines, see 
e.g. [US]. It is the goal of this paper to present a numer¬ 
ical approach to study the evolution of lipid membranes. 
We present several computations showing quite different 
shapes, and the influence of fluid flow on the membrane 
evolution. 

Since the classical papers of Canham [5] and Helfrich 
[6], there has been a lot of work with the aim of de¬ 
scribing equilibrium membrane shapes with the help of 
elastic membrane energies. Canham [5] and Helfrich [6] 
introduced a bending energy for a non-flat membrane, 
which is formulated with the help of the curvature of 
the membrane. In the class of fixed topologies the rele¬ 
vant energy density, in the simplest situation, is propor- 
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tional to the square of the mean curvature x. The re¬ 
sulting energy functional is called the Willmore energy. 
When computing equilibrium membrane shapes one has 
to take constraints into account. Lipid membranes have a 
very small compressibility, and hence can safely be mod¬ 
elled as locally incompressible. In addition, the presence 
of certain molecules in the surrounding fluid, for which 
the membrane is impermeable, leads to an osmotic pres¬ 
sure, which results in a constraint for the volume en¬ 
closed by the membrane. The minimal energetic model 
for lipid membranes consists of the Willmore mean cur¬ 
vature functional together with enclosed volume and sur¬ 
face area constraints. Already this simple model leads to 
quite different shapes including the biconcave red blood 
cell shapes, see [T]. 

Helfrich [6] introduced a variant of the Willmore en¬ 
ergy, with the aim of modelling a possible asymmetry of 
the bilayer membrane. Helfrich [6] studied the functional 
/(x — x)^, where x is a fixed constant, the so-called 
spontaneous curvature. It is argued that the origin of 
the spontaneous curvature is e.g. a different chemical en¬ 
vironment on both sides of the membrane. We refer to 
[7] and [8] for a recent discussion, and for experiments 
in situations which lead to spontaneous curvature effects 
due to the chemical structure of the bilayer. 

Typically there is yet another asymmetry in the bilayer 
leading to a signature in the membrane architecture. 
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This results from the fact that the two membrane layers 
have a different number of molecules. Since the exchange 
of molecules between the layers is difficult, an imbalance 
is conserved during a possible shape change. The total 
area difference between the two layers is proportional to 
M = j K. Several models have been proposed, which 
describe the difference in the total number of molecules 
in the two layers with the help of the integrated mean 
curvature. The bilayer coupling model, introduced by 
Svetina and coworkers inHn, assumes that the area per 
lipid molecule is fixed and assumes that there is no ex¬ 
change of molecules between the two layers. Hence the 
total areas of the two layers are fixed, and on assum¬ 
ing that the two layers are separated by a fixed distance, 
one obtains, to the order of this distance, that the area 
difference can be approximated by the integrated mean 
curvature, see ISHII]. We note that a spontaneous cur¬ 
vature contribution is irrelevant in the bilayer coupling 
model as this would only add a constant to the energy as 
the area and integrated mean curvature are fixed. 

Miao et al [12] noted that in the bilayer coupling model 
budding always occurs continuously which is inconsistent 
with experiments. They hence studied a model in which 
the area of the two layers are not fixed but can expand or 
compress under stress. Given a relaxed initial area dif¬ 
ference AHo, the total area difference AH, which is pro¬ 
portional to the integrated mean curvature, can deviate 
from AAq. However, the total energy now has a contri¬ 
bution that is proportional to (AH — AHq)^. This term 
describes the elastic area difference stretching energy, see 
Bin], and hence one has to pay a price energetically to 
deviate from the relaxed area difference. 

It is also possible to combine the area difference elas¬ 
ticity model (ADE-model) with a spontaneous curvature 
assumption, see Miao et al. [T2| and Seifert [T]. However, 
the resulting energetical model is equivalent to an area 
difference elasticity model with a modified AHq, see [T] 
for a more detailed discussion. 

It has been shown that the bilayer coupling model (BC- 
model) and the area difference elasticity model (ADE- 
model) lead to a multitude of shapes, which also have 
been observed in experiments with vesicles. Beside oth¬ 
ers, the familiar discocyte shapes (including the “shape” 
of a red blood cell), stomatocyte shapes, prolate shapes 
and pear-like shapes have been observed. In addition, the 
budding of membranes can be described, as well as more 
exotic shapes, like starfish vesicles. Moreover, higher 
genus shapes appear as global or local minima of the 
energies discussed above. We refer to mnMis] for more 
details on the possible shapes appearing, when minimiz¬ 
ing the energies in the ADE- and BC-models. 

Configurational changes of vesicles and membranes 
cannot be described by energetical considerations alone, 
but have to be modelled with the help of appropriate evo¬ 
lution laws. Several authors considered an -gradient 
flow dynamics of the curvature energies discussed above. 
Pure Willmore flow has been studied in USHU], where 
the last two papers use a phase field formulation of the 


Willmore problem. Some authors also took other aspects, 
such as constraints on volume and area miiiiis] , as 
well as a constraint on the integrated mean curvature 
naiio], into account. The effect of different lipid com¬ 
ponents in an L^-gradient flow approach of the curvature 
energy has been studied in [24H3Q] . 

The above mentioned works considered a global con¬ 
straint on the surface area. The membrane, however, 
is locally incompressible and hence a local constraint on 
the evolution of the membrane molecules should be taken 
into account. Several authors included the local inexten¬ 
sibility constraint by introducing an inhomogeneous La¬ 
grange multiplier for this constraint on the membrane. 
This approach has been used within the context of dif¬ 
ferent modelling and computational strategies such as the 
level set approach [311134] . the phase field approach [35L 
l37] . the immersed boundary method |38ll4Q] . the inter¬ 
facial spectral boundary element method [41] and the 
boundary integral method m- 

The physically most natural way to consider the local 
incompressibility constraint makes use of the fact that 
the membrane itself can be considered as an incompress¬ 
ible surface fluid. This implies that a surface Navier- 
Stokes system has to be solved on the membrane. The 
resulting set of equations has to take forces stemming 
from the surrounding fluid and from the membrane elas¬ 
ticity into account. In total, bulk Navier-Stokes equa¬ 
tions coupled to surface Navier-Stokes equations have to 
be solved. As the involved Reynolds numbers for vesicles 
are typically small one can often replace the full Navier- 
Stokes equations by the Stokes systems on the surface 
and in the bulk. The incompressibility condition in the 
bulk (Navier-) Stokes equations naturally leads to conser¬ 
vation of the volume enclosed by the membrane and the 
incompressibility condition on the surface leads a conser¬ 
vation of the membrane’s surface area. A model involving 
coupled bulk-surface (Navier-)Stokes equations has been 
proposed by Arroyo and DeSimone [43], and it is this 
model that we want to study numerically in this paper. 

Introducing forces resulting from membrane energies 
in fluid flow models has been studied numerically before 
by different authors, [3T14331 l36l l3711^ . However, typi¬ 
cally these authors studied simplified models, and either 
volume or surface constraints were enforced by Lagrange 
multipliers. In addition, either just the bulk or just the 
surface (Navier-)Stokes equations have been solved. The 
only work considering simultaneously bulk and surface 
Navier-Stokes equations are Arroyo et al [44] and Bar¬ 
rett et al. [45] [46] , where the former work is restricted 
to axisymmetric situations. In the present paper we are 
going to make use of the numerical method introduced 
in [46], see also [45] . 

The paper is organized as follows. In the next section 
we precisely state the mathematical model, consisting 
of the curvature elasticity model together with a cou¬ 


we introduce our numerical method which consists of an 
unfitted parametric finite element method for the mem- 


pled bulk-surface (Navier-)Stokes system. In Section HI 
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brane evolution. The curvature forcing is discretized and 
coupled to the Navier-Stokes system in a stable way us¬ 
ing the finite element method for the fluid unknowns. 
Numerical computations in Section IV demonstrate that 
we can deal with a variety of different membrane shapes 
and flow scenarios. In particular, we will study what in¬ 
fluence the membrane viscosity, the area difference elas¬ 
ticity (ADE) and the spontaneous curvature have on the 
evolution of bilayer membranes in shear flow. We finish 
with some conclusions. 


II. A CONTINUUM MODEL FOR FLUIDIC 
MEMBRANES 

We consider a continuum model for the evolution of 
biomembranes and vesicles, which consists of a curvature 
elasticity model for the membrane and the Navier-Stokes 
equations in the bulk and on the surface. The model is 
based on a paper by Arroyo and DeSimone [43] , where in 
addition we also allow the curvature energy model to be 
an area difference elasticity model. We first introduce the 
curvature elasticity model and then describe the coupling 
to the surface and bulk Navier-Stokes equations. 

The thickness of the lipid bilayer in a vesicle is typi¬ 
cally three to four orders of magnitude smaller than the 
typical size of the vesicle. Hence the membrane can be 
modelled as a two dimensional surface T in Given 
the principal curvatures xi and X 2 of T, one can define 
the mean curvature 


X = + X2 

and the GauB curvature 


Taking now into account that there is an optimal area 
difference AAq, the authors in [47H49| added a term pro¬ 
portional to 

(M(r) - Mof 

to the curvature energy, where Mq is a fixed constant 
which is proportional to the optimal area difference. 

For non-symmetric membranes a certain mean curva¬ 
ture Jc can be energetically favourable. Then the elastic¬ 
ity energy 0 is modified to 

/ (f ~ ds. 

The constant Jc is called spontaneous curvature. Taking 
into account that Jp ac K ds does not change for an 
evolution within a fixed topology class, the most general 
bending energy that we use in this paper is given by 
aE{T) with the dimensionless energy 

= 1 ( 2 ) 

where [3 has the dimension ( length )^- 

We now consider a continuum model for the fluid flow 
on the membrane and in the bulk, inside and outside 
of the membrane. We assume that the closed, time de¬ 
pendent membrane (r(t))t>o lies inside a spatial domain 
U C For all times the membrane separates Q into 
an inner domain U_(t) and an outer domain U+(t). De¬ 
noting by u the fluid velocity and by p the pressure, the 
bulk stress tensor is given by a = 2 pD^u) — p Id, with 
D{u) = I (V + (V u)'^) being the bulk rate-of-strain 
tensor. We assume that the Navier-Stokes system 

p (fit + (fl. V) fl) — V . a = 0 , V . fl = 0 


K = 

(as often in differential geometry we choose to take the 
sum of the principal curvatures as the mean curvature, 
instead of its mean value). The classical works of Gan- 
ham [5| and Helfrich |6] derive a local bending energy, 
with the help of an expansion in the curvature, and they 
obtain 

/ (f +aGK'j ds (1) 

as the total energy of a symmetric membrane. The pa¬ 
rameters a,aG have the dimension of energy and are 
called the bending rigidity a and the Gaussian bending 
rigidity ac • If we consider closed membranes with a fixed 
topology, the term Jp K ds is constant and hence we will 
neglect the Gaussian curvature term in what follows. 

As discussed above, the total area difference AA of the 
two lipid layers is, to first order, proportional to 

M(r) = J Kds. 


holds in 0_(t) and 0+(t). Here p and p are the density 
and dynamic viscosity of the fluid, which can take differ¬ 
ent (constant) values p±, p± in Arroyo and DeS¬ 

imone [43] used the theory of interfacial fluid dynamics, 
which goes back to Scriven [50] , to introduce a relaxation 
dynamics for fluidic membranes. In this model the fluid 
velocity is assumed to be continuous across the mem¬ 
brane, the membrane is moved in the normal direction 
with the normal velocity of the bulk fluid and, in addi¬ 
tion, the surface Navier-Stokes equations 

pr d* u - \/s. gr = iy E a fr , . fl = 0 

have to hold on T{t). Here pr is the surface material 
density, d* is the material derivative and . are the 

gradient and divergence operators on the surface. The 
surface stress tensor is given by 

= 2 Pr {u) — Pr 'Rr , 

where pr is the surface pressure, pr is the surface shear 
viscosity, Vr is the projection onto the tangent space and 

'^(u) = 5 2r (Vs M + (Vs up) 'Rr 




4 


is the surface rate-of-strain tensor. Furthermore, the 
term [a]l — a_ is the force exerted by the 

bulk on the membrane, where v denotes the exterior unit 
normal to The remaining term a fr denotes the 

forces stemming from the elastic bending energy. These 
forces are given by the first variation of the bending en¬ 
ergy aE{r{t)), see pi[43]. It turns out that fr points in 
the normal direction, i.e. fr = fr and we obtain, see 

mm, 

fr = -As K- {x-x) |Vs + ^{x-x)'^x 
+ ^(M(r)-Mo)(|Vsi7|2-x2) on T{t). 

Here is the surface Laplace operator, u is the Wein- 
garten map and |Vs + X 2 . Assuming e.g. no-slip 

boundary conditions on dft, the boundary of ^4, we ob¬ 
tain that the total energy can only decrease, i.e. 

<!»+<. Em) 

=-2 iJ.\Q{u)f dx + fir J |^s('w)|^ds^ <0. (3) 

We now non-dimensionalize the problem. We choose a 
time scale t, a length scale x and the resulting velocity 
scale u = x/i. Then we define the bulk and surface 
Reynolds numbers 


Re = X p-^u/jii-^ and Rer=xpru/pr^ 
the bulk and surface pressure scales 

p = P-^/t and pr = P-\-x/t = p-^u ^ 


and 


p* = p/p+ = 


.p-Ip+ 


in ^4+ 
in Q- 


p* = p/p+ = 

Pr = Pr/{p+ x), 


1 in 
A in ^4_ 


A = p_/p+, 


as well as the new independent variables x = xjx^t = tjt. 
For the unknowns 


a* = a/ (/i+ ux‘^) and x* = xJc^ Mq = Mq/x, /3* = x‘^ /3. 
We remark that the Reynolds numbers for the two re¬ 
gions in the bulk are given by Re and Rep*//i*, respec¬ 
tively, and that they will in general differ in the case of 
a viscosity contrast between the inner and outer fluid. 
In addition to the above equations, we of course also re¬ 
quire that u has zero divergence in the bulk and that the 
surface divergence of u vanishes on F. 

Typical values for the bulk dynamic viscosity p are 
around 10“^ —10“^^, see [4l|43l|52], whereas the surface 
shear viscosity typically is about 10“^ — 10“^^, see [H 
M US. The bending modulus a is typically 10 — 

10-19 see [301 EH ESI- 

The term = fXY/ (/li+ x) in Q suggests to choose the 
length scale 


x = ijlt/p+ Mr = l- 

As a* = aKfi^ux^) = x^) appears in Q, we 

choose the time scale 


i = p-^ x^/a . 


Choosing 


w = 5 . lO-"*? , , a = , 

s s m s^ 


see e.g. |43], we obtain the length scale 5 • 10“^m and the 
time scale 1.25s, which are typical scales in experiments. 
With these scales for length and time together with val¬ 
ues of ^ lO^kg/m^ for the bulk density and ^ 10“^kg/m^ 
for the surface densities, we obtain for the bulk and sur¬ 
face Reynolds numbers 


Re ^ 10 ^ and Rer ~ 10 ^ , 


and hence we will set the Reynolds numbers to zero in 
this paper. We note that it is straightforward to also 
consider positive Reynolds numbers in our numerical al¬ 
gorithm, see [45j [46] for details. Together with the other 
observations above, we then obtain the following reduced 
set of equations. 


u = u/u , p = p/p , pr = Pr/Pr 

we now obtain the following set of equations (on dropping 
the ^-notation for the new variables for ease of exposi¬ 
tion) 


—/i*A'u-}-Vp = 0 in^4±(t), 

- 2 Vs. ^s{u) + Vs. {pr 2r) 

= [2/i* — pid] ^ (a* ^ on r(t). (6) 


Re p* {ut E {u . u) — p* A u p = 0 in (t), 
Rer Pr d* u-Vs\^Pr ^s{u) - pr &) 

= [2 II* £(u) - pM] t + a* /r on r{t) , ( 4 ) 

with f* = f*u, 

/p = —X — (x — X*) |Vs d- px — x*)"^ X 

+ r (M(r) - Mo*) (IV on Tit) , (5) 


A downside of the scaling used to obtain is that the 
surface viscosity no longer appears as an independent 
parameter. However, studying the effect of the surface 
viscosity, e.g. on the tank treading to tumbling transition 
in shearing experiments, is one of the main focuses of this 
paper. It is for this reason that we also consider the fol¬ 
lowing alternative scaling, when suitable length and ve¬ 
locity scales are at hand. For example, we may choose the 
length scale x based on the (fixed) size of the membrane 
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and a velocity scale u based on appropriate boundary ve¬ 
locity values. In this case we obtain from Q, for small 
Reynolds numbers, the following set of equations 

— /i*A'u + Vp = 0 in (t), 

- Vs. (2 /ip 2s (u) - pr'Rr) 

= [2 D{u) — pl^~^ P ^ on r(t). (7) 


help of a weak formulation by Dziuk m, which is gener¬ 
alized by Barrett et al. [46] to take spontaneous curvature 
and area difference elasticity effects into account. A main 
ingredient of the numerical approach is the fact that one 
can use a weak formulation of that can be discretized 
in a stable way. In fact, defining A* = /3* (M(r) — Mg) 
and y = K-\- (A* — x*) the following identity, which has 
to hold for all y on F, characterizes /pi 


Note that here three non-dimensional parameters remain: 
/if, A and a*. Here /if compares the surface shear vis¬ 
cosity to the bulk shear viscosity, A is the bulk viscosity 
ratio and a* is an inverse capillary number, which de¬ 
scribes the ratio of characteristic membrane stresses to 
viscous stresses. Clearly, the system ^ corresponds to 
0 with /if = 1. Hence from now on, we will only con¬ 
sider the scaling 0 in detail. 

Of course, the system 0 needs to be supplemented 
with a boundary condition for u or a, and with an ini¬ 
tial condition for r(0). For the former we partition the 
boundary 90 of O into 9iO, where we prescribe a fixed 
velocity u = and ^ 20 , where we prescribe the stress- 
free condition an = 0, with n denoting the outer normal 
to O. 

We note that our non-dimensionalization may be dif¬ 
ferent to others presented in the literature. Often the 
length scale a = is chosen, where A(0) denotes 

the surface area of the vesicle at time zero, see e.g. 
[54] . Our length scale x may lead to simulations with 
A(0) = 47r5'^, with S > 0, so that our non-dimensional 
parameters in 0 correspond to the non-dimensional val¬ 
ues for a fixed length scale x = a as follows: 


M = 


Mr Mr — c'—* 

-= wr 5 a>c = b X 

b 


Ca = 


jjLj^ua'^ 


(a* 

( 8 ) 


III. NUMERICAL APPROXIMATION 

The numerical computations in this paper have been 
performed with a finite element approximation intro¬ 
duced by the authors in [46|. The approach discretizes 
the bulk and surface degrees of freedom independently. 
In particular, the surface mesh is not a restriction of the 
bulk mesh. The bulk degrees of freedoms u and p are 
discretized with the lowest order Taylor-Hood element, 
P2-P1, in our numerical computations. The evolution 
of the membrane is tracked with the help of paramet¬ 
ric meshes F^, which are updated by the fluid velocity. 
Since the membrane surface is locally incompressible, it 
turns out that the surface mesh has good mesh proper¬ 
ties during the evolution. This is in contrast to other 
fluid problems with interfaces in which the mesh often 
deteriorates during the evolution when updated with the 
fluid velocity, see e.g. [55] . 

The non-dimensionalized elastic forcing by the mem¬ 
brane curvature energy, 2 discretized with the 


(^fr^ x) = (Vs y, Vs x) + (Vs. y, Vs • x) 

- 2 ((Vs yf, £s(x) (Vs id)^) + (A* - X*) (x, [Vs x]^ *?) 

- i ([|x-x* - 2 (^. x)] Vs id, Vs x) 

- A* ((x. j7) Vs id, Vs x) • 

Here (•, •) is the L^-inner product on F, and VsV = 
with {ds-,,ds^,ds^Y = V^. Roughly speak¬ 
ing the above identity shows that /p has a divergence 
structure. We remark here that similar divergence struc¬ 
tures have been derived with the help of Noether’s theo¬ 
rem, see [saisTj. 

The numerical method of Barrett et al [46| has the 
feature that a semi-discrete, i.e. continuous in time and 
discrete in space, version of the method obeys a discrete 
analog of the energy inequality 0. In addition, this 
semi-discrete version has the property that the volume 
enclosed by the vesicle and the membrane’s surface area 
are conserved exactly. After discretization in time these 
properties are approximately fulfilled to a high accuracy, 
see Section |TV| The fully discrete system is linear and 
fully coupled in the unknowns. The overall system is 
reduced by a Schur complement approach to obtain a 
reduced system in just velocity and pressure unknowns. 
For this resulting linear system well-known solution tech¬ 
niques for finite element discretizations for the standard 
Navier-Stokes equations can be used, see Barrett et al. 

m- 


IV. NUMERICAL COMPUTATIONS 

In shearing experiments the inclination angle of the 
vesicle in the shear flow direction is often of interest. Here 
we will always consider shear flow in the xi direction with 
^3 being the flow gradient direction. Precisely, if Q = 
[— L, I/]^ X [—VF, IF], then we prescribe the inhomogeneous 
Dirichlet boundary condition g{x) = (x3,0,0)^ on the 
top and bottom boundaries dift = [—x {±W}. 
Assuming the vesicle’s centre of mass is at the origin, then 
— ~ /o (t) ld — x0x dx denote the vesicle’s moment 

of inertia tensor. Let p, with |^ = 1 and pi > 0, be the 
eigenvector corresponding to the smallest eigenvalue of 
M. Then the vesicle’s inclination angle is defined by ^ = 
arg(pi + ip 3 ) G (—7r/2,7r/2], where arg : C ^ (—7r,7r]. 
For later use we also note that the deformation parameter 
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D is defined by (6 — c)/(6 + c), where 6, c are the major 

and minor semiaxes of an ellipsoid with the same moment 

1 

of inertia tensor, see e.g. [59]. Hence, in 2d, D = (Amax — 
111 

^min)/(-^max + where Amax and Amin are the two 

eigenvalues of M. 

The inclination angle 0 is important for the classifi¬ 
cation of different types of dynamics in the shear flow 
experiments that we will present. The classical defor¬ 
mation dynamics for vesicles are the tank treading (TT) 
and the tumbling (TU) motions. In the tank treading 
motion the vesicle adopts a constant inclination angle 
in the flow, while the surface fluid rotates on the mem¬ 
brane surface. This motion is observed for small vis¬ 
cosity contrasts between the inner and the outer fluid 
and, as we will see later, at low surface membrane vis¬ 
cosity. At large viscosity contrasts or large membrane 
viscosity the tumbling motion occurs. In the tumbling 
regime the membrane rotates as a whole, and the incli¬ 
nation angle oscillates in the whole interval (—7r/2,7r/2]. 
In the last ten years a new dynamic regime for vesicles 
in shear flow has been identified. In this regime the in¬ 
clination angle is neither constant nor does it oscillate 
in the whole interval (—7r/2,7r/2]. The dynamics are 
characterized by periodic oscillations of the inclination 
angle 0 such that 0 G [—Oq^Oq] for a Oq in the open in¬ 
terval (0,7r/2). This regime was first predicted theoret¬ 
ically by m and subsequently observed experimentally 
in m- Later this regime has been studied by differ¬ 
ent groups, see e.g. (53] (54] l62ll66] for more details. In 
[60] this motion was called vacillating-breathing and later 
the same motion was also called trembling, transition 
mode or swinging. Following [64] we will refer to this 
new regime as the transition (TR) mode. 

In our numerical simulations we will only consider the 
scaling 0. For all the presented simulations we will 
state the reduced volume as a characteristic invariant. It 
is defined as = 6 7r^ V(0)/^®(0), see e.g. [15]. Here 
V{t) and A{t) denote the volume of the discrete inner 
phase and the discrete surface area, respectively, at time 
t. Moreover, if nothing else is specified, then our numer¬ 
ical simulations are for no-slip boundary conditions, i.e. 
diCl = dCl and ^ = 0. In all our experiments it holds 
that X* /3* = 0, and we will only report the values of x* 
and for simulations where they are nonzero. Here we 
recall, as stated in the introduction, that the energy 

for X* /3* 7 ^ 0 is equivalent to with X* = 0, the same 
value of P* > 0, and a modified value of Mq . Finally, we 
stress that our sign convention for curvature is such that 
spheres have negative mean curvature. 

A. 2d validation 

In order to validate our numerical method, we repro¬ 
duce some numerical results from naEzi, where we al- 



FIG. 1: (Color online) A plot of Oj it against Ar for 
L = 20, kF = 5, A = 1, Of* = 0.01, Re = 10“^; compare 
with [32] Fig. 1]. 



FIG. 2: (Color online) A plot of D against A for 
L = 20, kF = 5, A = 1, Of* = 0.01, Re = 10“^; compare 
with [32] Fig. 2]. 


ways consider a domain ft = [—L, L] x [—kF, kF]. As these 
works consider Navier-Stokes flow in the bulk, we con¬ 
sider Q with Re = 10“^, Rep = 0, //p = x* = /3* = 0 
and vary A. For the comparison with Figures 1-3 in [32] 
we also set a* = 10“^. Moreover, we consider vesicles 
with reduced areas Ar = ^ {0.6,0.7,0.8,0.9}, 

and with a = = 1, so that perimeter and area are 

given by P(0) = 2 ir and Fl(0) = Ar^r. At first, for 
A = 1, we try to recreate m Fig. 1]. To this end, we 
set L = 20 and kF = 5, and use stress-free boundary 
conditions left and right, rather than periodic bound¬ 
ary conditions in the -direction on the square domain 
[—5, 5]^ as used in [32]. We obtain the results in Figure 
where we plot O/tt against Mr, which show a good agree¬ 
ment with [32] Fig. 1]. Similarly, in trying to recreate 
[32] Fig. 2] we also compute the deformation parame¬ 
ter D, and plot D against the excess length parameter 

A = 2 (1 — Ar)/{7T A?). We obtain the results in Fig¬ 
ure which show good agreement with [32] Fig. 2]. In 
Figure [^ we plot the critical viscosity ratio Ac for the TT 
to TU transition against the reduced area Mr- It should 
be noted that our numerical method produces larger val- 
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FIG. 3: (Color online) A plot of Ac against Ar for 
L = 20, IF = 5, Of* = 0.01, Re = 10“^; compare with 
m Fig. 3]. 


FIG. 5: (Color online) Analogue of the phase diagram 
from [54l Fig. 8] for the domain ft = [—3,3]^, starting 
with a prolate shape with Vr = 0.8 and the longest axis 
in the xs direction. 



FIG. 4: (Color online) A plot of ^/f, V/^ and 
against A for Ar = 0.8, L = 11.55, IF = 3.85, a* = 2; 
compare with [611 Fig. 1]. 


B. 3d validation 

For a similar validation in 3d we compare our method 
to some numerical results from m, where Stokes flow in 
an infinite domain is considered. In order to reproduce 
the phase diagram in [54l Fig. 8], which also contains nu¬ 
merical results from [66], we let ft = [—3,3]^ and choose 
the initial shape of the interface to be a prolate vesicle 
with a reduced volume of Vr = 0.8 and a surface area of 
^(0) = 477, so that S = 1. The results from our algo¬ 
rithm are shown in Figure Due to finite size effects, 
and the different boundary conditions, we observe differ¬ 
ent critical values for the phase transitions compared to 
[54l Fig. 8]. However, qualitatively our numerical method 
produces similar results. 


ues of Ac than reported in |32j Fig. 3]. 

Moreover, in trying to recreate EH Fig. 1], we also ran 
with Re = 0.05, L = 11.55 and IF = 3.85, so that the 
restriction parameter y as defined in [67] is y = 0.26. 
However, we note that periodic boundary conditions in 
the -direction are used in m, with the length L of the 
domain not clearly stated. We obtain the results in Fig¬ 
ure]^ where apart from 0, normalized by f, we also show 
the membrane tank treading velocity V = -p^ /p \u\ ds, 
normalized by and the tumbling frequency cj, normal¬ 
ized by ^ (note that the frequency in [67| Fig. 1] is said 
to be normalized by ^). It should be noted that qual¬ 
itatively our results agree well with [67] Fig. 1], but our 
numerical method produces a smaller value of Ac than 
reported in im Fig. 1], 

Overall we are satisfied that our numerical method per¬ 
forms well. The observed differences with existing results 
in the literature can be explained by differences in the 
length of the domain, different boundary conditions and 
different numerical methods used. 


C. Effect of surface viscosity 

We consider the effect that surface viscosity has on the 
TT to TU transition. To this end, we let ft = [—3,3]^, 
and choose as initial shape of the vesicle a biconcave 
shape with reduced volume Vr = 0.8 and M(0) = 477, 
so that S = 1. We let A = 1. In Figure [6a| we present a 
phase diagram with the axes labelled in terms of the non- 
dimensional values Ga = ^ and Ai = recall ® for 
S = 1. The evolutions for a* = 0.1, and either /ip = 3, 
/ip = 1, or /ip = 0.1, are visualized in Figure [61^ where we 
observe the motions TU, TR and TT, respectively. We 
stress that the tumbling occurs for a viscosity contrast 
of A = 1, and so is only due to the chosen high surface 
viscosity /ip. The fact that vesicles undergo a transition 
from steady tank treading to unsteady tumbling motion 
has been observed earlier by [68], where, however, the 
authors used a particle-based mesoscopic model to ana¬ 
lyze the fluid vesicle dynamics. A plot of the inclination 
angle 0 for the simulations in Figure [^ can be seen in 
Figure [^ 

























0 







(a) Phase diagram. (b) Top to bottom: — 3, 1, 0.1 for 1/a* — 10 at times t = 0, 3, 6, 9. 


FIG. 6: (Color online) Phase diagram for A = 1 for the domain Vt = [—3,3]^, starting with a biconcave shape with 
Vr = 0.8 and the shortest axis in the xi-direction. The three big circles in the phase diagram correspond to the 

simulations in (b). 



FIG. 7: The inclination angle 0 for the computations in 
Figure [61^ They correspond to the motions TU, TR 
and TT, respectively. 



(a) K = —5 (b) X = 5 



FIG. 9: (Color online) Phase diagram for 
Ca = 1/a* = 10 for the domain ft = [—3,3]^, starting 
with biconcave shapes with Vr = 0.8 and the shortest 
axis in the -direction. 


FIG. 8: (Color online) The vesicles for x* = ±5 at time 

t = 0. 


D. Effect of spontaneous curvature 


Here the initial shapes of the vesicles, for a reduced 
volume of Vr = 0.8 and surface area ^(0) = 47r, so that 
S' = 1, were chosen to be numerical approximations of lo¬ 
cal minimizers for the curvature energy Jp(x — x*)^ ds. 
These discrete local minimizers were obtained with the 
help of the gradient flow scheme from m, and for the 
choices x* = ±5 they are displayed in FigureFor Ca = 
l/of* = 10 we show a phase diagram of = /if versus 
X* in Figure where the initial vesicles are aligned such 
that their shortest axis is in the -direction. Similarly, 
in Figure we show a phase diagram of = /if ver¬ 
sus X* when the initial vesicles are aligned such that 
their shortest axis is in the X 2 -direction. The results in 
Figures and IT indicate that the values of the surface 



FIG. 10: Phase diagram for Ca = l/a* = 10 for the 
domain ft = [—3,3]^, starting with biconcave shapes 
with Vr = 0.8 and the shortest axis in the X 2 -direction. 
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/xf 

X* = —5 

X* = 0 

X* = 5 

0.05 

0.179 

0.178 

0.194 

0.1 

0.169 

0.158 

0.179 

0.2 

0.161 

0.116 

0.138 


TABLE I: Some inclination angles 0 for the TT motions 
in Figure]^ 


* 

Mr 

1 

II 

o 

II 

II 

0.05 

0.180 

0.156 

0.162 

0.1 

0.160 

0.132 

0.156 

0.2 

0.118 

0.084 

0.143 


TABLE II: Some inclination angles 0 for the TT 
motions in Figure pT} 


small time step sizes. 

In our next simulation, we let Vt = [—2.5, 2.5]^ and set 
A = Of* = 1, as well as /3* = 0.46 and Mg = —33.5. As 
initial vesicle we take a varying-diameter cigar-like shape 
that has Vr = 0.75 and ^(0) = 9.65, so that S = 0.88. A 
simulation can be seen in FigureAs a comparison, we 
show the simulation with = 0 in Figure Similarly 
to previous studies, where an energy involving area dif¬ 
ference elasticity terms was minimized, we also observe 
in our hydrodynamic model that less symmetric shapes 
occur when the ADE-energy contributions are taken into 
account. 


G. Shearing for budded shape (two arms) 


viscosity, at which the transitions between TT, TR and 
TU take place, strongly depend on the spontaneous cur¬ 
vature as well as on the orientation of the initial vesicle. 
As in [68], where the case x* = 0 was studied, we also 
observe that the inclination angle in the tank treading 
motion decreases as /if increases, see Tables [I| and [11} 


E. Shearing for a torus 


Here we use as the initial shape a Clifford torus, that is 
aligned with the X 2 -x^ plane, with reduced volume Vr = 
0.71 and ^(0) = 13.88, so that S = 1.05._We let A = 
/if = 1, a* = 0.05 and use the domain ft = [—2,2]^. 
See Figure El where the torus appears to tumble whilst 
undergoing strong deformations. 

Repeating the experiment for an initial torus aligned 
with the shear flow direction, and setting a* = 1 and 


/if = 10, leads to the results shown in Figure 12 This 
shows a TR motion. Setting /if =0, on the other hand, 
leads to TT, as shown in Figure [l^ A plot of the incli¬ 
nation angle 0 for the simulations in Figures and 13 
can be seen in Figure [T^ while we visualize the flow in 
the xi-xs plane in Figure [Ts] 


F. Effect of area difference elasticity 

We consider O = [—4,4]^ and set A = /if = a* = 1. 
The parameters for are /3* = 0.053 and Mg = —48.24. 
For the vesicle we use a cup-like stomatocyte initial shape 
with Vr = 0.65 and .4(0) = 82.31, so that S = 2.56. See 
Figure [T6| for a numerical simulation. As a comparison, 
we show the same simulation with /d* = 0 in Figure [TT] 
In Figure we show the evolutions of the discrete 
volume of the inner phase and the discrete surface area 
over time. Clearly these two quantities are preserved al¬ 
most exactly for our numerical scheme in this simulation. 
In fact, a semidiscrete variant of our scheme conserves 
these two quantities exactly, and so in practice the fully 
discrete algorithm will preserve them well for sufficiently 


We start a scaled variant of the final shape from Fig¬ 
ure in a shear flow experiment inQ = [—2, 2]^. In par¬ 
ticular, the initial shape is axisymmetric, with reduced 
volume Vr = 0.75 and .4(0) = 5.43, so that S = 0.66. 
We set A = /if = 1, Of* = 0.05. See Figure for 
a run with = 0.1 and Mg = —33.5. We observe 
that the shape of the vesicle changes drastically, with 
part of the surface growing inwards. This is similar to 
the shapes observed in Figure [^ where the presence 
of a lower reduced volume led to cup-like stomatocyte 
shapes. We repeat the same experiment for = 0 in 
Figure [^ Now the budding shape loses its strong non¬ 
convexity completely, as can be clearly seen in the plots 


Plots of the 
where 


23 


of the two-dimensional cuts in Figure pT 
bending energy a* F^*(r^) are shown in Figure 
we recall that the energy inequality in ^ does not hold 
for the inhomogeneous boundary conditions employed in 
the present simulations. 


H. Shearing for a seven-arm starfish 

We consider simulations for a scaled version of the final 
shape from Barrett et al [H Fig. 23] with reduced vol¬ 
ume Vr = 0.38 and .4(0) = 10.54, so that S = 0.92, inside 
the domain Q = [—2, 2]^. We set A = /if = a* = 1. In or¬ 
der to maintain the seven-arm shape during the evolution 
we set /3* = 0.05 and Mg = 180. The first experiment 
is for no-slip boundary conditions on dQ and shows that 
the seven arms grow slightly, see Figure [^ If we use 
the shear flow boundary conditions, on the other hand, 
we observe the behaviour in Figure [^ where we have 
changed the value of a* to 0.05. The vesicle can be seen 
tumbling, with a tumbling period of about 7, with the 
seven arms remaining intact throughout. Repeating the 
same experiment with ^* = 0 yields the simulation in 
Figure!^ Not surprisingly, some of the arms of the vesi¬ 
cle are disappearing. We also tried to investigate whether 
the arms enhance or inhibit the tumbling behaviour of 
the vesicle. To this end, we repeated the simulation in 
Figure!^ for an ellipsoidal vesicle with the same reduced 
volume and surface area. This vesicle also exhibited TU 
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FIG. 11: (Color online) Shear flow for a torus with A = /ip = 1. The plots show the interface F^ within 11, as well as 
cuts through the xi-xs plane, at times t = 0, 2.5, 5, 7.5. The interface at t = 10 is very close to the plot at t = 2.5. 



FIG. 12: (Color online) Shear flow for a torus with A = 1, /ip = 10. The plots show the interface F^ within H, as 
well as cuts through the xi-x^ plane, at times t = 0, 2.5, 5, 7.5. 


with a tumbling period of about 7, so there was no sig¬ 
nificant change to the behaviour in Figure 

V. CONCLUSIONS 

We have introduced a parametric finite element 
method for the evolution of bilayer membranes by cou¬ 
pling a general curvature elasticity model for the mem¬ 
brane to (Navier-) Stokes systems in the two bulk phases 
and to a surface (Navier-)Stokes system. The model 
is based on work by Arroyo and DeSimone [43], which 
we generalized such that area difference elasticity effects 
(ADE) are taken into account. Our main purpose was to 
study the influence of the area difference elasticity and of 
the spontaneous curvature on the evolution of the mem¬ 
brane. In contrast to most other works, we discretized 
the full bulk (Navier-) Stokes systems coupled to the sur¬ 
face (Navier-) Stokes system and for the first time coupled 
this to a bending energy involving ADE and spontaneous 
curvature. 

The numerical simulations led to the following findings. 

• The proposed numerical method conserves the vol¬ 
ume enclosed by the membrane and the surface area 
of the membrane to a high precision. 

• The transition from a tank treading (TT) motion to 
a transition motion (TR) and to a tumbling (TU) 
motion depended strongly on the surface viscosity. 


We observed that the surface viscosity alone with 
no viscosity contrast between inner and outer fluid 
can lead to a transition from tank treading to a 
TR-motion and to tumbling. Similar observations 
have been reported by [68] using a particle-based 
method. 

• The surface viscosity at which a transition be¬ 
tween the different motions TT, TR and TU oc¬ 
cur, strongly depends on the spontaneous curvature 
and on the initial alignment of the vesicle. In par¬ 
ticular, we observed that for negative spontaneous 
curvature and an initial biconcave vesicle aligned 
such that the shortest axis is in the shear flow di¬ 
rection all transitions occurred for larger values of 
the surface viscosity. Eor this alignment, and for 
positive spontaneous curvature, we observed that 
tumbling occurred already for much smaller values 
of the surface viscosity. The reverse was true for an 
alternative alignment. Here we recall that our sign 
convention for curvature means that spheres have 
negative mean curvature. 

• In some cases, shear flow can lead to drastic shape 
changes, in particular for the ADE-model. Eor ex¬ 
ample, we observed the transition of a budded pear¬ 
like shape to a cup-like stomatocyte shape in shear 
flow if an ADE-model was used for the curvature 
elasticity. 

• The ADE-model can also lead to starfish-type 
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FIG. 13: (Color online) Shear flow for a torus with A = 1, /if = 0. The plots show the interface F^ within 11, as well 

as cuts through the xi-xs plane, at times t = 0, 2.5, 5, 7.5. 



FIG. 14: The inclination angle 0 for the simulations in 
Figures 


12 


and 


13 



shapes with several arms, see e.g. [DUS]. In com¬ 
putations for a seven-arm starfish for a model in¬ 
volving an ADE type energy, we observed that in 
shear flow the overall structure seems to be quite 
robust. In particular, the seven arms deformed but 
remained present even in a tumbling motion. How¬ 
ever, arms tend to disappear if the area difference 
elasticity term is neglected. 

Thus we have shown that the proposed numerical method 
is a robust tool to simulate bilayer membranes for quite 
general models which in particular take the full hydrody¬ 
namics and a curvature model involving area difference 
elasticity and spontaneous curvature into account. 


FIG. 15: (Color online) The flow at time t = 7.5 in the 
xi-xs plane for the simulations in Figures 


12 


and 


13 
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FIG. 16: (Color online) Flow for a cup-like stomatocyte shape with Vr = 0.65 for Mg = —48.24 and /3* = 0.053. The 
plots show the interface F^ at times t = 0, 5, 10, 20, with the top row visualizing the triangulations by showing half 

the vesicle. 




FIG. 18: The evolutions of the relative discrete volume 
V(t)/V(0), and the relative discrete surface area 
A{t)/A{^) over time. 
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FIG. 19: (Color online) Flow for a varying-diameter cigar-like shape with Vr = 0.75 for Mg = —33.5 and f3* = 0.46. 
The plots show the interface F^ at times t = 0, 1, 10, 50, with the top row visualizing the triangulations. 




FIG. 20: (Color online) Same 


as Figure 19 with /3* = 0. 


University Regensburg, Germany. 

[47] U. Seifert, L. Miao, H.-G. Dobereiner, and M. Wortis, 
in The Structure and Conformation of Amphiphilic Mem¬ 
branes, Springer Proceedings in Physics, Vol. 66, edited 
by R. Lipowsky, D. Richter, and K. Kremer (Springer- 
Verlag, Berlin, 1992) pp. 93-96. 

[48] W. Wiese, W. Harbich, and W. Helfrich, J. Phys. Con- 
dens. Matter 4, 1647 (1992). 

[49] B. Bozic, S. Svetina, B. Zeks, and R. Waugh, Biophys. 
J. 61, 963 (1992) 

[50] L. E. Scriven, Chem. Eng. Sci. 12, 98 (1960) 

[51] T. J. Willmore, An. §ti. Univ. “Al. I. Cuza” Ia§i SecU I 
a Mat. (N. S.) IIB, 493 (1965). 

[52] M. Eaivre, Drops, vesicles and red blood cells: Deforma- 
bility and behavior under flow, Ph.D. thesis, Universite 
Joseph Eourier, Grenoble, Erance (2006). 

[53] D. Abreu, M. Levant, V. Steinberg, and U. Seifert, Adv. 
Colloid Interface Sci. 208, 129 (2014) 

[54] A. P. Spann, H. Zhao, and E. S. G. Shaqfeh, Phys. Eluids 
26, 031902 (2014) 

[55] E. Bansch, Numer. Math. 88, 203 (2001) 

[56] R. Capovilla and J. Guven, J. Phys. Condens. Matter 16, 
2187 (2004) 


[57] D. Lengeler, (2015), arXiv:1506.08991 

[58] J. W. Barrett, H. Garcke, and R. Niirnberg, J. Sci. 
Comp. 63, 78 (2015) 

[59] S. Ramanujan and C. Pozrikidis, J. Eluid Mech. 361, 117 
(1998) 

[60] C. Misbah, Phys. Rev. Lett. 96, 028104 (2006) 

[61] V. Kantsler and V. Steinberg, Phys. Rev. Lett. 96, 
036001 (2006) 

[62] N. J. Zabusky, E. Segre, J. Deschamps, V. Kantsler, and 
V. Steinberg, Phys. Eluids 23, 041905 (2011) 

[63] H. Noguchi and G. Gompper, Phys. Rev. Lett. 98, 128103 
(2007) 

[64] A. Yazdani and P. Bagchi, Phys. Rev. E 85, 056308 

( 2012 ) 

[65] A. Earutin and C. Misbah, Phys. Rev. Lett. 109, 248106 

( 2012 ) 

[66] T. Biben, A. Earutin, and C. Misbah, Phys. Rev. E 83, 
031921 (2011) 

[67] B. Kaoui, T. Kruger, and J. Harting, Soft Matter 8, 
9246 (2012) 

[68] H. Noguchi and G. Gompper, Phys. Rev. E 72, 011901 
(2005) 






































14 



FIG. 21: (Color online) Shear flow for a budding shape with A = /rf = 1. Here /3* = 0.1 and Mg = —33.5. The plots 
show the interface F^ within 11, as well as cuts through the xi-xs plane, at times 
t = 0, 5, 15, 17.5, 20, 25, 27.5, 32.5. 



FIG. 22: (Color online) Same as Figure 21 but with f3* = 0. 
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FIG. 23: The bending energy a* £’*(r^) for the 
computations in Figures 


21 


and 


22 




FIG. 24: (Color online) Flow for a seven-arm figure 
with Vr = 0.38. Here /3* = 0.05 and Mq = 180. The 
triangulations F^ at times t = 0 and 5. 



FIG. 25: (Color online) Shear flow for a budding shape with A = /if = 1. Here /3* = 0.05 and Mq = 180. The plots 
show the interface F^ within Q at times t = 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5. 
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FIG. 26: (Color online) Same as Figure 


with = 0. 
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